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We consider the distribution of the (properly normalized) 
numbers of nodal domains of wave functions in 2-d quantum 
billiards. We show that these distributions distinguish clearly 
between systems with integrable (separable) or chaotic un- 
derlying classical dynamics, and for each case the limiting 
distribution is universal (system independent). Thus, a new 
criterion for quantum chaos is provided by the statistics of 
the wave functions, which complements the well established 
criterion based on spectral statistics. 
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The attempts to discern the "finger-prints" of the un- 
derlying classical dynamics in the patterns observed in 
the quantum wave-functions is one of the main themes 
of quantum chaos 10-^]. Here, we present a new ap- 
proach which is based on the distribution of the number 
of the nodal domains of wave- functions. We show that 
this distribution makes a clear distinction between inte- 
grable (separable) and chaotic systems, and we study the 
universal features associated with each of the cases. As a 
convenient paradigm we consider "quantum billiards" in 
2-d whose classical dynamics correspond to geodesic flows 
with specular reflections on the domain boundaries. The 
corresponding Schrodinger (Hclmholtz) equation reads: 



A*(r) = EV(r) ; r e ft 



(1) 



where ft is a connected compact domain on a 2-d Rieman- 
nian manifold, and A is the corresponding Laplacian. If 
ft has boundaries, Dirichlet boundary conditions are as- 
sumed throughout. 

Considering only real solutions of ([j]), the nodal lines 
are the zero sets of ^(r). The nodal domains are the 
connected domains where ^(r) has a constant sign. The 
distribution of the number of nodal domains is the sub- 
ject of the present note. Courant was the first to 
address this counting problem. Ordering the eigenfunc- 
tions by the magnitude of the corresponding eigenvalues, 
and denoting by v n the number of nodal domains of the 
n'th eigenfunction, Courant proved that 



v n < n 



(2) 



Pleijel H showed that the equality in (||) holds only for 
a finite sequence of eigenf unctions, and that 



-ImT^ <(- 
Ji 



0.691 , 



(3) 



Here, we associate to each wave-function the normal- 
ized nodal-domain number 



n 



< £„ < 1 



(4) 



and study their distribution, defined in the following way. 
Consider the spectral interval I g (E) = [E,E+gE], (g > 
arbitrary) , and denote by Nr the number of eigenvalues 
in I g (E). The £ distribution associated with I g {E) is 

E n ei g (E) 

We shall show the existence of the limiting distribution, 



P(0= lim P(Z,I g (E)) 



(6) 



and present its universal features. 

For domains with boundaries, we also study the dis- 
tribution of the number of nodal intersections with the 
boundary, denoted by v n . This is the number of times 
the normal derivative at the boundary vanishes. In this 
case, the appropriate normalized parameter is 



fin 



(7) 



where j\ is the first zero of the Bessel function Jq(x). 



and we define the distribution P(r), I g (E)) and the lim- 
iting distribution P(rj) as in (j^) and (^), respectively. 

We shall first discuss integrable systems, and in par- 
ticular, those generic integrable systems which are also 
separable This set of systems includes the rectan- 
gular and elliptic domains in R 2 , surfaces of revolution 
Q and Liouville surfaces In the two latter cases, 

we shall deal only with surfaces for which exist global 
action-angle variables for the geodesic flows. The rele- 
vant common features to the systems mentioned above 
are that (a) The hamiltonian H{I\,I<2) is a homoge- 
neous function of degree 2 of the action variables I\ % 
(which are expressed in units of H). (b) The contour 
line r defined by H(Ii,l2) = 1 in the positive quad- 
rant of the (iijia) plane is monotonic. In the semiclas- 
sical limit of interest here, the spectrum is given by the 
EBK quantization: E^ m = H(l + ai,m + a 2 ) + 0(\f~E) 
where I, m are the integer quantum numbers and ai,2 
are the Maslov indices. The normalized number of nodal 
domains (Q) is £/. m = vi im /m im . Seperability implies 
vi,m = I m + 0(1), and = J\[(Ei >m ), where Af(E) is 
the spectral counting function. To leading order, AT(E) 
can be replaced by the first term in the Weyl series, 
JV(E) = AE(1 + 0(^7=)), where A is the area enclosed 



1 



by r and the I\, I2 axes. Substituting these approximate 
values in (^|), and converting the discrete sum into an 
integral, incure a relative error of order 0(^=). Thus, 



1 



(8) 
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dhdi 2 s e 



H(i 1 ,i 2 )ei g 



h h 



AH(h,I 2 ) 
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The homogeneity of the hamiltonian function leads to an 
expression for P(£, I g (E)) which is independent of both 
E and g. This establishes the existence of the limit- 
ing distribution (||). Changing the integration variables 
(7i,/ 2 ) — ► (£,s), where £ — P(ii,/ 2 ) and s is the arc- 
length along the line T, (0) is reduced to 



ds S £ 



A 



(9) 



This form allows a simple geometrical interpretation: 
Ii ( s )l2 (s) is the area of the rectangle generated by a point 
(ii, I2) on r and its projections on the axes. Thus, P(f) 
is the probability distribution of the areas of rectangles 
scaled by A, the area between T and the axes. How- 
ever, r is monotonic, and therefore the scaled areas are 
all smaller than 1. Thus, P(£) = for £ > £ m , where 
£ m is the maximum value of the scaled areas. The mono- 
tonicity of T ensures the existence of an explicit function 
for T: I2 = 12(^1)1 hi terms of which the integral (||) can 
be computed for every value of £ < £ m . The resulting 
distribution is 



p(0 







1 v J 2 (Ii)-4(7i)Ji 

2 / 2 (/i)+^(/i)/i 
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(10) 



The sum is over all the real values of /1 which satisfy 
£ _ J i J 2( 7 i) j n ^]^c vicinity of £ m , typically two solutions 
coalesce, leading to a square root singularity of P(£) near 
£ m . The vanishing of P(£) in the interval (£ m , 1), and 
the square root singularity at £ m are the universal fea- 
tures which characterize the nodal domain distributions 
for separable, integrable systems. 

The distribution P(rj) can be derived in a similar way, 
the only difference being that v n is twice the value of one 
of the action variables I\ , say. (The rectangular billiard is 
the only exception, where v n = 2(ii + I2), for I1.2 > 1). 
The scaling with respect to \fn is naturaly called for. 
Denoting by l[ the intersection of T with the 1% axis, 

the maximum value of 77 is rj m = 2/-^ j\fA. Then, the 
analogue of ( |Tol) reads, 







V < Vrn 



(HI 



Near r? m , P(rf) reflects the way that T intersects the I\ 
axis, and it is not universal. 




FIG. 1. The nodal domains distributions (|9|) for the rect- 
angular and circular billiards in the spectral intervals: rectan- 
gle- 62500 < n < 125000; circle- 30000 < n < 60000. Smooth 
line- the limiting distributions (lid). 



Numerical simulations for the rectangular and the cir- 
cular domains are shown in Figs. (|l|, ^]), together with 
the predicted limiting distributions. The deviations can 
be explained by considering the next to leading contribu- 
tions in \j\f~E. The singularities seen in the numerical 
data mark the contributions of periodic orbits |1 1| . 




FIG. 2. The distribution of the normalized numbers of 
boundary intersections for the rectangular and circular bil- 
liards, in the same spectral intervals as in Fig. 1. Smooth 
line- the limiting distributions (|ll|). 



The most important feature of wavefunctions corre- 
sponding to separable systems is that the nodal set forms 
a simple grid of intersecting nodal lines. In chaotic sys- 
tems, nodal intersections may occur, and when q nodal 
lines intersect at a point, they do so at angles The 
intersections are rare, and are unstable against small per- 
turbations. There is no known analytical expression for 
the number of nodal domains, and therefore we cannot 
offer an expression for the limiting distribution of the nor- 
malized numbers of nodal domains for chaotic billiards. 
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However we shall present physical arguments augmented 
by numerical evidence which support the conjecture that 
the nodal domains of random waves reproduce the ob- 
served distributions of chaotic billiards M. While this 



work was in progress, Bogomolny and Schmit |16| pro- 
posed a method to compute the nodal domain distribu- 
tions for random waves, based on a percolation model. 
We shall compare our numerical result to their predic- 
tions. 

The scattering approach to the quantization of billiards 
p2| , p^[ implies that the wavefunction can be well approx- 
imated by 

L 

¥„(r) « ai J '( fc « r ) eM > a - 1 = "H- 1 )'- ( 12 ) 

l=-L 

Here, k = VE, L = [k£/2Tt] and C is the boundary 
length. The 2L + 1 coefficients {a;} are the eigenvector 
of the scattering matrix S(k — k n ) corresponding to an 
eigenvalue 1. The S matrix is the semiclassical quanti- 
zation of the billiard bounce map. For hyperbolic bil- 
liard maps, quantum ergodicity implies that the vectors 
{a{\ are uniformly distributed on the sphere. For large 
L, this distribution can be approximated by considering 
{Re(a;), Im(a;)} as independent random Gaussian vari- 
ables with (|a;| 2 ) = 1. Hence, the set of eigenfunctions 
( |l2| ) are conjectured to be simulated by the Gaussian 
random wave ensemble. 

To check the conjectured relation between wave func- 
tions of chaotic billiards and the ensemble of random 
waves, we computed the lowest 1637 (1483) eigenfunc- 
tions of the stadium (Sinai) billiards. The nodal do- 
mains of the billiard wave functions were counted using 
the Hoshen-Kopelman method [119]. We performed the 
following tests. 

(i.) The distribution of boundary intersections. Given 
a function of the type ( |l2| ) we computed the number of 
the zeros of its normal derivative on a circle of radius R 
(L = kR). Denoting the normal derivative by u(9) = 
d r ^(r, 9)\ r= R the number of its zeros is fll4| 



d(9 



d£dr; 



(1 



) (13) 
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(14) 



of nodal intersections limits to a delta function centered 
at T]^^- = 1. The agreement between the billiards data 
shown in Fig. (^|) and the predictions of the random 
waves model, lends very strong support to our conjec- 
ture. However, the fluctuations in v n are affected by the 
presence of "bouncing ball" states, which explains the 
apparent difference between the widths of the 77 distribu- 
tions for the two billiards. 



where u(9) = dgu(6). The mean number of zeros (v), and 
its variance were evaluated by performing the Gaussian 
integrals, with the result 



The expression for (i>) is exact, while that for the variance 
is the leading term in the semiclassical limit. The vari- 
ance was computed analytically up to an integral whose 
numerical value is quoted. Since for a billiard of area A, 

n — 4fr7 ( PI ) supports the proposed scaling (Q). More- 
over, it implies that the distribution of the scaled number 




FIG. 3. The number of boundary intersections u n as a func- 
tion of n for the stadium and the Sinai billiards. Smooth line 



- the random waves result p4[). Inset - the n distribution. 

(ii.) Signed area distribution. Denote by |f2(^)|± the 
total area where the wave function is positive (negative) . 
Clearly (|f2|±) = |^|/2, where |f2| is the billiard area. To 
compute the signed area variance, ((|f2|+ — |0|_) )/|£l| 2 , 
we used the identity 



|fi(*)l± = ^- Km / d 2 r [ d£ 



(15) 
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FIG. 4. The normalized signed area variance for the sta- 
dium and Sinai billiards. Smooth curve - analytic expression 
for random waves. Triangles - numerical simulation for ran- 
dom waves (see text). 

For a circle of radius R the signed area variance can 
be reduced to an expression involving a single integral 
which was evaluateded numerically, and is shown in Fig. 
4. as a function of n = (Rk/2) 2 . For large n, 
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0.0386n" 



(16) 



Fig. 4. shows also the data for the stadium and the 
Sinai billiards, which converge to the asymptotic limit 
( |l6| ) . The slow convergence of the Sinai data is due to the 
exsitence of corners with sharp angles in the desymetrized 
Sinai billiard. The effect of such corners on the wave 
function is accenuated at low energies, and we checked 
this effect, by computing the variance in a random wave 
model for which the only angular waves allowed are mul- 
tiples of 4 (triangular dots). This simulates the boundary 
conditions on the symmetry line of the Sinai billiard, and 
shows similar slow convergence to the asymptotic value. 



1000- 



Sinai (xlOO) 




ory does not include the effects of the boundaries, and 
therefore it applies in the limit of high n. It is consistent 
with our numerical data, and confirms the linear scaling 
with n. The numerical simulations based on the random 
waves model, supported by the results of [[u| predict that 
var(v n ) scales linearily with n for large n. Thus, the lim- 
iting distribution for chaotic billiards is expected to 
limit to a delta function, centered at the value predicted 
by |1(| to be 0.0624, much less than the Pleijel bound, 
or the values observed for the square and circle billiards. 

To summarize, we provided here ample evidence to 
show that the study of nodal domains counts provides 
a clear criterion of quantum chaos, which is not directly 
related to spectral statistics. Our results add a new di- 
mension to previous attempts to attribute the morphol- 
ogy of wave functions to the nature of the underlying 
classical dynamics |l| || . 
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FIG. 5. The number of nodal domains v n (solid lines) and 
inner nodal domains (dotted) for the Sinai and the stadium 
billiards, and for random waves (circles for nodal domains 
and squares for interior nodal domains). Smooth lines: the 
prediction of the theory of Bogomolny and Schmit. 

(Hi.) The number of nodal domains v n . The number of 
nodal domains for the stadium and Sinai billiards are 
plotted in Fig. (5). The predicted (v n ) for the ran- 
dom waves model are shown as large circular dots. These 
data were computed numerically by sampling the random 
wave ensemble, assuming a square domain. The agree- 
ment between the billiards and the random waves results 
is satisfactory, and confirms the expected scaling v n ~ n 
for high n. The random and the billiards data sets, how- 
ever, show clear deviations from the expected asymptotic 
scaling for low n values, where v n ps n 1 / 2 . To check the 
reason for this behaviour we computed and plotted the 
number of interior nodal domains - those domains whose 
boundaries do not include the billiard boundaries. Their 
number for the billiards and the random waves model 
scale as n, suggesting the following explanation: At low 
n, most of the domains are boundary domains whose 
number scales as the number of boundary intersections, 
that is, n 1 / 2 . Only when the wavelengths are sufficiently 
smaller than the typical linear dimension of the billiard, 
the majority of domains are interior, and their number 
increases linearily with n. The smooth lines in Fig 5. 
represent the analytical result of |l(| for (v n ). This the- 
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